source ../../TMD.tcl
source entropy.tcl
# interface 
proc calcEnthalpy { temp w1 w2 } {
   set nw1 [llength $w1]
   set nw2 [llength $w2]
   set c_w1 [doublesToByteArray $w1]
   set c_w2 [doublesToByteArray $w2]
   return [c_enthalpy $temp $nw1 $c_w1 $nw2 $c_w2]
}
proc ave { list } {
   set n [llength $list]
   set c_list [doublesToByteArray $list]
   return [c_stat $n $c_list]
}
proc min { list } {
  set min [lindex $list 0]
  foreach v $list {
    if { $min > $v} { set min $v}
  }
  return $min
}
# c code 
critcl::ccode {
    #include "math.h"
    #include <stdlib.h>
    #include "stdio.h"
    }
critcl::cproc testd { double n} double {
 printf("exp(%f) = %e\n",n, exp(n));
 return exp(n);
}
critcl::cproc c_enthalpy { double temp int nstate1 double* state1 int nstate2 double* state2 } double {
   double energy=0.0,weight=0.0, beta;
   int i,j,k;
   double w1,w2,w_ave,ws;
   double* d= (double*) malloc(nstate1*nstate2*sizeof(double));
   double* s= (double*) malloc(nstate1*nstate2*sizeof(double));
   beta=1.98e-3*temp;
   k=0;
   for (i=0;i< nstate1; i++) {
      for (j=0; j< nstate2; j++) {
        w1=*(state1+i);
        w2=*(state2+j);
        *(d+k)=w1-w2;
        s[k]=w1+w2;
        //printf("%f %f %f %f\n",w1,w2,d[k],s[k]);
        k++;
      }
   }
   w_ave=0.0;
  for (i=0;i< k;i++) {
   w_ave+=s[i];
  }
  w_ave/=k;
  for (i=0;i< k;i++) {
   ws=s[i]-w_ave;
  if (fabs(ws) < 500.0) {
    
     //printf("%f,%f\n",d[i],ws);
     energy+=d[i]*exp(-beta*ws);
     weight+=exp(-beta*ws);
  }
  }
  energy/=weight;
  //printf("%e\n",w_ave);
  free(d);
  free(s);
  return energy;
}
critcl::cproc c_stat { int num double* list} double {
  double ave=0.0,std=0.0;
  double min=list[0],max=list[0];
  int i;
  for ( i=0; i<num ;i++) {
    ave+=list[i];
    std+=list[i]*list[i];
    if (list[i]< min) min=list[i];
    if (list[i]> max) max=list[i];
  }
  ave/=num;
  std/=num;
  std-=ave*ave;
  return ave;
}
#list<->bytearray
 proc listToByteArray { valuetype list {elemsize 0} } {
    if { $valuetype == "i" || $valuetype == "I" } {
       if { $::tcl_platform(byteOrder) == "littleEndian" } {
          set valuetype "i"
       } else {
          set valuetype "I"
       }
    }

    switch -- $valuetype {
    f - d - i - I {
       set result [binary format ${valuetype}* $list]
    }
    s {
       set result {}
       foreach elem $list {
          append result [binary format a$elemsize $elem]
       }
    }
    default {
       error "Unknown value type: $valuetype"
    }
    }

    return $result
 }

 interp alias {} stringsToByteArray {} listToByteArray s
 interp alias {} intsToByteArray    {} listToByteArray i
 interp alias {} floatsToByteArray  {} listToByteArray f
 interp alias {} doublesToByteArray {} listToByteArray d

 #
 # Generic routine to convert a bytearray into a list
 #
 proc byteArrayToList { valuetype bytearray {elemsize 0} } {
    if { $valuetype == "i" || $valuetype == "I" } {
       if { $::tcl_platform(byteOrder) == "littleEndian" } {
          set valuetype "i"
       } else {
          set valuetype "I"
       }
    }

    switch -- $valuetype {
    f - d - i - I {
       binary scan $bytearray ${valuetype}* result
    }
    s {
       set result  {}
       set length  [string length $bytearray]
       set noelems [expr {$length/$elemsize}]
       for { set i 0 } { $i < $noelems } { incr i } {
          set elem    [string range $bytearray \
                         [expr {$i*$elemsize}] [expr {($i+1)*$elemsize-1}]]
          set posnull [string first "\000" $elem]
          if { $posnull != -1 } {
             set elem [string range $elem 0 [expr {$posnull-1}]]
          }
          lappend result $elem
       }
    }
    default {
       error "Unknown value type: $valuetype"
    }
    }
    return $result
 }

 interp alias {} byteArrayToStrings {} byteArrayToList s
 interp alias {} byteArrayToInts    {} byteArrayToList i
 interp alias {} byteArrayToFloats  {} byteArrayToList f
 interp alias {} byteArrayToDoubles {} byteArrayToList d


# data structure 
proc select { vlist n } {
    set level(1) $vlist
    for {set l 2} { $l <= $n } {incr l } {
	set level($l) {}
	foreach group $level([expr $l -1]) {
	    set idx [expr 1+[lsearch $vlist [lindex $group end]]]
	    foreach item [lrange $vlist $idx end] {
		if {[lsearch $group $item] < 0} {
		    #puts "$group: $item"
		    lappend level($l) [concat $group $item]
		}
	    }
	}
    }
    return $level($n)
}
proc flaten { list } {
    return $list
}
# time
proc time { block {ntime 1}} {
    set t1 [clock clicks -milliseconds]
    for { set i 0} { $i < $ntime } { incr i } {
	set result [ uplevel 1 $block]
    }
    set t2 [clock clicks -milliseconds]
    set second [expr double($t2-$t1)/(1000*$ntime)]
    #return $second
    puts "Execute: $block in $second seconds\n Result:\n $result"
}



